{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression modeling with statsmodel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load example data\n",
    "spector_data = sm.datasets.spector.load()\n",
    "spector_data.exog = sm.add_constant(spector_data.exog, prepend=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.416\n",
      "Model:                            OLS   Adj. R-squared:                  0.353\n",
      "Method:                 Least Squares   F-statistic:                     6.646\n",
      "Date:                Thu, 23 Jun 2016   Prob (F-statistic):            0.00157\n",
      "Time:                        17:57:48   Log-Likelihood:                -12.978\n",
      "No. Observations:                  32   AIC:                             33.96\n",
      "Df Residuals:                      28   BIC:                             39.82\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.4639      0.162      2.864      0.008         0.132     0.796\n",
      "x2             0.0105      0.019      0.539      0.594        -0.029     0.050\n",
      "x3             0.3786      0.139      2.720      0.011         0.093     0.664\n",
      "const         -1.4980      0.524     -2.859      0.008        -2.571    -0.425\n",
      "==============================================================================\n",
      "Omnibus:                        0.176   Durbin-Watson:                   2.346\n",
      "Prob(Omnibus):                  0.916   Jarque-Bera (JB):                0.167\n",
      "Skew:                           0.141   Prob(JB):                        0.920\n",
      "Kurtosis:                       2.786   Cond. No.                         176.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# fit and summarize OLS model\n",
    "mod = sm.OLS(spector_data.endog, spector_data.exog)\n",
    "res = mod.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More detailed example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create artificial data\n",
    "np.random.seed(9876789)\n",
    "nsample = 100\n",
    "x = np.linspace(0, 10, 100)\n",
    "X = np.column_stack((x, x**2))\n",
    "beta = np.array([1, 0.1, 10])\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we need an intercept so create a column of ones\n",
    "X = sm.add_constant(X)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 4.020e+06\n",
      "Date:                Thu, 23 Jun 2016   Prob (F-statistic):          2.83e-239\n",
      "Time:                        17:57:48   Log-Likelihood:                -146.51\n",
      "No. Observations:                 100   AIC:                             299.0\n",
      "Df Residuals:                      97   BIC:                             306.8\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.3423      0.313      4.292      0.000         0.722     1.963\n",
      "x1            -0.0402      0.145     -0.278      0.781        -0.327     0.247\n",
      "x2            10.0103      0.014    715.745      0.000         9.982    10.038\n",
      "==============================================================================\n",
      "Omnibus:                        2.042   Durbin-Watson:                   2.274\n",
      "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.875\n",
      "Skew:                           0.234   Prob(JB):                        0.392\n",
      "Kurtosis:                       2.519   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# fit and summary\n",
    "model = sm.OLS(y, X)\n",
    "results = model.fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [  1.34233516  -0.04024948  10.01025357]\n",
      "R2:  0.999987936503\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', results.params)\n",
    "print('R2: ', results.rsquared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# now simulate artificial data with a non-linear relationship between x and y\n",
    "nsample = 50\n",
    "sig = 0.5\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))\n",
    "beta = [0.5, 0.5, -0.02, 5.]\n",
    "\n",
    "y_true = np.dot(X, beta)\n",
    "y = y_true + sig * np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.933\n",
      "Model:                            OLS   Adj. R-squared:                  0.928\n",
      "Method:                 Least Squares   F-statistic:                     211.8\n",
      "Date:                Thu, 23 Jun 2016   Prob (F-statistic):           6.30e-27\n",
      "Time:                        17:57:48   Log-Likelihood:                -34.438\n",
      "No. Observations:                  50   AIC:                             76.88\n",
      "Df Residuals:                      46   BIC:                             84.52\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.4687      0.026     17.751      0.000         0.416     0.522\n",
      "x2             0.4836      0.104      4.659      0.000         0.275     0.693\n",
      "x3            -0.0174      0.002     -7.507      0.000        -0.022    -0.013\n",
      "const          5.2058      0.171     30.405      0.000         4.861     5.550\n",
      "==============================================================================\n",
      "Omnibus:                        0.655   Durbin-Watson:                   2.896\n",
      "Prob(Omnibus):                  0.721   Jarque-Bera (JB):                0.360\n",
      "Skew:                           0.207   Prob(JB):                        0.835\n",
      "Kurtosis:                       3.026   Cond. No.                         221.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# fit and summary\n",
    "res = sm.OLS(y, X).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [ 0.46872448  0.48360119 -0.01740479  5.20584496]\n",
      "Standard errors:  [ 0.02640602  0.10380518  0.00231847  0.17121765]\n",
      "Predicted values:  [  4.77072516   5.22213464   5.63620761   5.98658823   6.25643234\n",
      "   6.44117491   6.54928009   6.60085051   6.62432454   6.6518039\n",
      "   6.71377946   6.83412169   7.02615877   7.29048685   7.61487206\n",
      "   7.97626054   8.34456611   8.68761335   8.97642389   9.18997755\n",
      "   9.31866582   9.36587056   9.34740836   9.28893189   9.22171529\n",
      "   9.17751587   9.1833565    9.25708583   9.40444579   9.61812821\n",
      "   9.87897556  10.15912843  10.42660281  10.65054491  10.8063004\n",
      "  10.87946503  10.86825119  10.78378163  10.64826203  10.49133265\n",
      "  10.34519853  10.23933827  10.19566084  10.22490593  10.32487947\n",
      "  10.48081414  10.66779556  10.85485568  11.01006072  11.10575781]\n"
     ]
    }
   ],
   "source": [
    "print('Parameters: ', res.params)\n",
    "print('Standard errors: ', res.bse)\n",
    "print('Predicted values: ', res.predict())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFrCAYAAADilOGNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xlc1OX2wPHPIOQKrrmholKmiEuLppUGmWkuZVkuoV0s\n67ZomXnb7kU0W6/aZvXrVpamZovZYqVpGomprWqCWi6ICu7bDCiyzPn98YgCojLD7Jz36zUvcvzO\n93kGjDPPdo5FRFBKKaWU9wR5uwNKKaVURafBWCmllPIyDcZKKaWUl2kwVkoppbxMg7FSSinlZRqM\nlVJKKS8LdncDFotFz04ppZSqUETE4sj1bg/GAHqW2X9NmDCBCRMmeLsbygn6s/Nv+vPzXxaLQ3EY\n0GlqpZRSyus0GCullFJepsFYnVNMTIy3u6CcpD87/6Y/v4rF4u71XIvFIrpmrJRSqqKwWCy+uYGr\nNM2bNyc9Pd1bzfu0iIgItm/f7u1uKKWU8hCvjYxPfnJwa9v+Sr83Sinlv5wZGeuasVJKKeVlGoyV\nUkopL9NgrJRSSnmZBuMyGjFiBOPHj/d2N5RSSgUgDcYuFhsby3vvveftbiillPIjGoyVUkopL/O5\nYJyWls6wYROJjU1k2LCJpKU5fhbZFfdYs2YNl19+OTVr1mTIkCHk5OQAcOTIEfr370/9+vWpW7cu\n/fv3JzMzE4D//Oc/JCcnM2rUKMLCwnjooYcAGDNmDM2aNaNmzZp06tSJFStWONwfpZRSAUxE3Pow\nTZyptOe3bdsukZGPCmQJiECWREY+Ktu2bS/1HqVxxT1yc3MlIiJCXn31VcnPz5d58+ZJSEiIJCQk\nyKFDh2T+/PmSk5MjWVlZMmjQIBkwYMCp18bExMj06dOL3W/OnDly+PBhKSgokJdeekkaNmwoJ06c\nOGv7Z/ueKaWU8n0nf4c7FisdfYHDDTgQjOPiJhQJonIqmMbFTSjzN8EV91i+fLmEh4cXe+6qq66S\nhISEM65ds2aN1KlT59SfSwvGJdWuXVv+/PPPs/69BmOllPJfzgRjn5qmzsiwA9VLPFudzEy7R++R\nmZlJeHh4seciIiIAOH78OP/85z9p3rw5tWrV4tprr+XIkSPnzJg1ZcoUoqKiqF27NrVr18ZqtXLg\nwIEy90cppVRg86lgHB4eBGSXeDabxo3L3k1X3KNRo0ZkZGQUe27Hjh0ATJ06lc2bN/Prr79y5MgR\nli9fDnAqGJcsKr1ixQomT57MvHnzOHz4MIcPHyYsLEzTXSqllDrFp4LxpEnxREYmcjqYZhMZmcik\nSfEevUfXrl0JDg5m2rRp5OfnM3/+fH755RcAbDYbVatWJSwsjEOHDjFhwoRir23QoAHbtm079Web\nzUZISAh169YlNzeXp59+GpvNVua+KKWUCnw+FYxbtIhgyZLRxMVNITY2kbi4KSxZMpoWLSI8eo+Q\nkBDmz5/P+++/T926dfn0008ZOHAgAI888gjHjh2jXr16XHXVVfTp06fYax9++GE+/fRT6taty5gx\nY+jduze9evWiVatWtGjRgmrVqtG0adMy90UppVTg06pNPki/N0op5b+0apNSSinlhzQYK6WUUl6m\nwVgppZTyMg3GSimllJdpMFZKKaXKKzcX5s6Fbt2cerkGY6WUUspZJ07Af/4DzZrBu+/CI484dZtg\nF3dLKaWUqjguuACqVIEffoA2bZy+zXnPGVsslulAP2CviLQ/+dx/gf7ACWArMEJErGd5vZ4zdpB+\nb5RSyn+565zx+0CvEs8tBtqKSEdgM/CkI40qpZRSfmXRInjzTbfd/rzBWERWAIdLPPe9iBSWQVoN\nNHFD37wiNDSUsLAwwsLCqFSpEtWqVTv13Ny5c73dPaWUUp60dy8MHQoPPggXXeS2ZlyxgesuYKEL\n7uMTbDYbVqsVq9VKREQE33zzzannhg4desb1BQUFXuilUkopt7Lb4Z13oF07iIiA9evhhhvc1ly5\nNnBZLJZ/A3ki8uG5rita2SgmJoaYmJjyNOsxhUWfi0pISGDz5s0EBQXxzTffMG3aNJYsWcLFF1/M\n+PHjAVi6dCkjR44kLS0NgIyMDEaPHs2KFSsIDQ3l0Ucf5YEHHvD4+1FKKVVG//43LFsGS5ZAhw7n\nvDQpKYmkpKRyNed0MLZYLPFAH+C6811bssygv/viiy+YP38+H374ITk5OSxZsuSMawrrGosI/fr1\nY8iQIXz66aekp6dz/fXX06ZNG2JjYz3ddaWUUmXx+OPwzDNQqdJ5Ly05yJw4caLDzZU1GFtOPswf\nLJbewL+A7iJywuFW/dw111xzqnRilSpVznntypUrsdlsPP744wC0bNmSu+66i48++kiDsVI+JC0t\nnYSEGWRk2AkPD2LSpHiHSq8iAqmp0KqVOe5S0rRp5hd7rVpQu7b5GhUFNWu67D0oF6pVy6PNnTcY\nWyyWD4EYoK7FYtkBJAJPARcAS06OAFeLiEvnXS0ObQo/O3ecEHKkHvGOHTtIT0+nTp06J/sj2O12\nDcRK+ZC0tHR69pzG1q0TgepANqtXJ56/FvqePfD997B4sflaubL574svPvPaI0dg9244fNj896FD\nsHkzbN8OYWFuemfqvPLzwWqFk7+jveW8wVhE7ijl6ffd0JcS7bq7BedZSnxSqF69OseOHTv15927\nd5/676ZNm9KqVStSU1M91j+llGMSEmYUCcQA1dm6dSIJCVOYPTux9BeNGgVz5kBsLPTsCePHQ2Tk\n2UcSCQlnPpedDdWrn/m88owdOyAuDq6+Gl54watd0XSYLtCxY0e++eYbjhw5wu7du5k2bdqpv+va\ntSsXXHABL730EidOnKCgoICUlBT++OMPL/ZYKVVURoad04G4UHUyM+2lXW6MHw/798P8+XD//ebY\ni6NTemcLxMuXw913m81D+fmO3VOVzbx5cMUV0K8fPPect3ujwfhcSo6AzyY+Pp7WrVsTERFBnz59\nih2BqlSpEt9++y2//PILzZs3p379+tx3333YbDZ3dVsp5aDw8CAgu8Sz2TRuHHT2abr69SHYTRmF\nL7kEoqPhqaegdWv48ENz1EaVX3Y23Huv2aD19dfma5D3Q+F502GWuwFNh+kw/d4o5VmlrRlHRiay\n4sm2NHx5qlkPbtjQO51btswcs8nKgs8+MxvElPPeeANWrzZf3bRW70w6TA3GPki/N0p5XuFu6sxM\nOx1Cj/BM9gaq70g3u6B7lcwI7GEi8O23cN11ULWqd/vi70Rct0P4LDQYBwj93ijlJcePm408r78O\n48bB2LFmh7RSDnAmGGsJRaWUKrR7tzlutGaNqU/rD3780Zxbbt/e2z3xPXv3QoMG3u5FmXh/1Vop\npXxFy5Zms5S/BGKAjAxztOq++8zubmXODY8YAQMG+PY52SI0GCullD+74w7YtMmsJUdFwUsvQW6u\nt3vlPcnJJpd0SIg5Gubm9WFX0TVjH6TfG6U8wG43v6j95Jd1mWzaBI88Yo5cLVjg7d54Vm4uJCbC\njBnwv//BTTe59PaOpEvVNWOllCqLY8dg+HAYONCMLANF69awcCEcPOjtnnje2rXmw8i6deYMuAs5\nnS7VAToy9kH6vVHKjfbsMaOmSy6Bd9/V3dLqvIYNm8icOeMonqUtm7i40tOlOjMy1jVjpVTFkZIC\nXbqYFIgffFCxAnFWFkyeDJr9z2FOpUt1kAbjc5gxYwbt27enevXqNG7cmAceeICjR48Cpl7l8OHD\nS33dihUruPrqq6lVqxb16tWjW7du/P77757sulKqpB9/NEkznnvO5JUOpLXissjONlO5F11kgnJ2\nyfSffuKPP+CVVzza5DnTpbqIBuOzmDp1Kk8++SRTp07FarWyevVq0tPTueGGG8g/mbi9tNzVNpuN\n/v378/DDD3P48GEyMjJITEykckX6BK6UL2re3BR1CKQ1Ykc0aGCqTC1bBr/8YoLyK6+YRCf+IDMT\n4uOhb1+oUcOjTU+aFE9kZCKnA7JJlzppUrzL2tA141LYbDYaN27MjBkzGDhw4Knns7OzadmyJS++\n+CLp6els3bqVDz74oNhrf//9d3r27MmhQ4ecbt+XvzdKqQCxbh1MmAAPPwwxMd7uzdkdOwZTpsCr\nr8I995jiGV6o/1w0XWrjxq7fTe2bwdhmM2s70dEQGupcw+W4x3fffUf//v3JyckhqEQ1j/j4ePLy\n8mjVqlWpwdhms9GyZUv69u3LkCFD6NKlC7Vq1XKofQ3GSil10uOPQ1oavPgitGjh7d6USWBs4LLZ\noFs36N7dfHVms0E573HgwAHq1at3RiAGaNSoEQcOHDjra0NDQ1mxYgVBQUHce++91K9fn5tvvpn9\nmhlHKeUvdu6ETz4xo1Jve+450xc/CcTO8r1gnJICqammoPa6dWY6ovBgvsViplVKM2HC6WvCwsxr\n8/NhwwZzPwfUq1ePAwcOYC+lfuju3bupV6/eOV9/ySWX8N5777Fjxw5SUlLIzMxkzJgxDvVBKVUO\n8+fDs896uxf+6/BhmD4dwsPhzjvN2eW8PPe0JWI2ZU2ZUvrfV6rknnZ9jO8F4+hoaNvWpDLr0MHk\nGBU5/ThXMC68xmo9nQ4tKsrczwFdu3alcuXKzJ8/v9jzWVlZLFy4kB49epT5Xq1atSI+Pp6UlBSH\n+qCUctLnn8P998ONN3q7J/6rfXv47juTRKNzZ5g0yQTmL75wzf3z82HLFjP1HB0Nt91mfm8XFLjm\n/n7I9zJwhYaa3KKpqSaIOrNmXM57hIWFMX78eEaPHk1oaCg9evRg165dPPjggzRr1ozhw4fz3HPP\nUVBQwIkTJ069zmKxkJaWxjfffMPgwYMJDw9n586dzJ07l65duzr+PpRSjvn8c1MwYeFCuOwyb/fG\n/zVoAKNGmce2bVClSunXvfceHDlippJbtDCB+8ABaNq09J3PXbua5Ct9+pjUlVdfHRBHzQo3eTlF\nRNz6ME2c6WzP+5L33ntPoqOjpVq1atKwYUO5//775ciRIyIiMmHCBAkKCjr1sFgs0rRpU8nMzJRB\ngwZJeHi41KhRQ5o0aSL333+/2Gy2MrfrD98bpXzO55+L1K8v8ttv3u5JxTN3rshDD4n07y8SHS1S\nu7ZIq1Yiv/9e+vUFBZ7tnwds27ZdIiMfFcgq/B3uUKz0zd3UFZx+b5RyUH4+xMbCyy/DFVd4uzeq\nAiqeMlMLRSilKqLgYFi+PCCmOpV/Kj1lZtlpMFZKBQYfCcSOlNpTgaMwZWYN7GQ58XqdpvZB+r1R\nyj+VVmovMtK1pfaUb0pLS2dAj6nMTEviUtYHQNIPpZQ6n61bzTFGH5OQMKNIIAaoztatE53fYav8\nRosgSL7qMB1Y79TrNRgrpfzLqlXmaMymTd7uyRk8UWpP+ZikJJPp8bLLCAsJxhIZ6dRtdM1YKeU/\nfv0Vbr7Z1CJu08bbvTnD6VJ7xYvQu7LUnvKi0moe1KkDjz0GvXrBBReYa5woZKFrxj5IvzdKlWLN\nGujdG959F/r393ZvSqVrxgHs0CFTD7swmVRy8lkTSvlV1abmzZuTnp7u1rb9VUREBNu3b/d2N5Ty\nHampcP318PrrUKSsqS9ypNSe8gMFBSbD2L//bXJ25+ebVMvLl0OXLqW+xK+CsVJKldnevfDLLz47\nIlYB6ocfYMwYM+387LPw0EOm+FBUVLGRcW4urFgBX39tHps3azBWSimvOXDADJh+XWaj3p4UcltF\nExYeSt26UK8ep77WqwdV811Qt125x65dMHo0rF0Lkyeb2RiLxawHn5ym3p8TysKFJvguWQKtWkG/\nfuZx2WWagUsppTxm71748cfTj507oUdnG/+X0o16+1LZd2Fb/tsvmeMbv2eDzcK+rGrssVXn6BHh\n/wruIVK2YG3Slkork6nVVAOyVxXdnFWpElx5JcydW6w4xp7sUN7/oQsLxp5eOenXD157DRo2LF/z\nOjJWSqkyEoFvP7aR8lEKn2yIZtv+UG66IpPBdb+ns3UJtaY9Q/C+TOjevfja4nvvwb59kJ0Nx46R\nm5FBSHo6FiCXYHpX/ZHaN15FXBz07QuVK3v7nVYwNps5nnSWzVmHD5sB8v/+B7ffDrfeCtdee/af\nkzNrxrrfXinlWxYuhCef9HYvzrBqFVx/pY2L7urGuAXdWb2/JYcatGbmH+3oU7CAejdfQ3CdsOI1\n2Qvrqb/9tqkFvGQJabM/pFPQTayjHScIYQNt2Hfh53TqdJBp06B5oxM8NnQnSUlgP2ozDdtspfYp\nLS2dYcMmEhubyLBhE0lLq1ibYl32/lNSTCDOzzdrwqmpgPns9PzzZgp63z4za/3WW3DDDW74wORo\nmSdHH2g5QKVUWX30kSmDuHKlt3tySlqayODBIk2aiHybsFLswcEiIBIUJDJ9ukh+/pkvslpFVq0y\nX0uIi5sgkCU1sMqVrJIaWAWyJC5ugoiI7F60VrKr1ZWNVTrI3qAGkh8ULPntOpxxr+Il+0QgSyIj\nH5Vt27a749sg27Ztl7i4CRITM17i4ia4rR1H+lPu979ggcjtt4scPSrSoYNISIhIhw5y4oBVXn9d\npFEjkUGDRDZtcqxvOFFCUYOxUso3vP22SOPGIuvWebsnImJ+Pz/+uEidOiITJ4pkZYkJiEV+aZcW\nbM8nJmb8yeBR/BEbO/70RXl5Iq+/LnZLkAhIPkGy7MFPi5UBLgzqxe9zOqi7kqcDf1mU6/0fPSpy\n990iLVqI/Pijec5qlfwVq2TOW1Zp3lykd++zl2M+H2eCsU5TK6W8b/JkeO45k1qwffsyv8wd07T5\n+WYqslUrCNr8F+v/FMaPh+rVIe3AIUZe0ocHoocz8pI+pB045PD9T2fpKqpElq7gYLjzTizt20FI\nCAV1LmTxYguXXw7ff28ucWXqzfN9H30x57bT73/5cujQweyOXrfOrO8DS1aH0uGfXXhzVigzZ5rV\nkssuc0vXS+do9Hb0gY6MlVLncuyYmSrcudOhl7ljtPbddyJRUSI3X7VP9t12n8iFF4r8/bdL23Po\nPkWmu+12kU8/FWnZUuTGG0X69n3DJSPjsvSnTKN5D3NoZGy1mqWPhQvN3POCBaf+6tAhkfh4kWbN\nRL78UsRuL3/f0GlqpVRF4cpp2oICkYQEkfbhB2RLv4fEXqeOyEMPiRw86Jb2CtdfY2MdX389cULk\nlVdE6tbNl8Y1vpfLSJYurJQa7Hbqw0FZ3pcnp8TLqswfagqXFoKDRdq3N5sATvrsMxObH3zQqRWH\ns3ImGOs5Y6WUX3LVNG1WFtx5J1yQ9hdr93fA8k2umaN+5plix1tcOS3cokUEs2cnOvw6MLUIHn4Y\n7ryzEjPjjvLgwkEEYSczrB4Fn/9AcwdTb5blfU2aFM/q1Yln5NyeNGm0U+/BFVq0iGDJktEkJEwp\nknq0lBzgRXdKb9wIe/awp0pzRo2C9evh44/NqSZv02CslPIsEbNeV06uqJCUng433QSXXw5vTd2H\n5YZ8079t28wv8CK5h32tIlPt2jAmoRGyxIIlX2hi3c/e8dPgw5ehatUy36cs76vMgc/Dzvmh5sQJ\nc/6o8KjZhg1Imyg+XNuWR26Cu++G2bOL5fTwLkeH0o4+0GlqpZSIWYx75RUzJ+gC5V3DTU4WadhQ\n5KWXTq4TnmentC/uKC7ssz0kRI40bi0LL+gvB0KbizX9UJlv4ZPvqzwKCkSmTRNp1crM6YuIWK2S\nOX+VDOhhlY4dnd8lXVY4MU2tGbiUUu539CjcdZcZin7yCbRs6ZLbOlshafp0k1dk1ixThvaUIrmH\nS8sX7ZMVmYr0+VBeKK/c9Scz/mjPm2+aVI1l4ZPvyxk7d5p/ZzabqXndqhVbt6YzZEgqa9deS9u2\nq/nkk4to1cq9702rNimlfM+aNSaHYK9eMHWqV+cF8/Nh3Dj45as9fDCnEhd1vdBrfXGnZcvg3nvN\n9Purr0LD6gFclMJmM4u/qammzOHDD8Pjj0NwMEuXZtCv3wlyciKASniqvrSmw1TKT1SYNIZJSSZ3\n4DPPwBtveDUQHz4MffpA3urfWZHbmYu2LPJaX9ztuutMfGrRArpG2zgY1Q3p1g2uvvqsqTX9UmFO\n6WuvNeUN58+Hf/+bfIJ58UXo27c2OTnhmEAMvnA++mx0A5dSHpaWlk7PntOK7Uxdvdr9n9a94sor\nYeVKuPhir3Zj5wYb43qn8I/Wf3HH1n9heestUxYvgFWtCi+8AHe1SSEsPhULBcj69VgmT4bERFOZ\nyN8V3SltsUBwMOvWmZnqOnWgY8f/8fPPj5R4kXM74N1NR8ZKeZgvZjNym6pVvR6Id6TasHW8hrm7\nriEu6V4sX34Z8IG4qFa3RhPcvi0FlUJIC4pk+zuLKWjfEb75xuzX8mdFinLYW0cxaX5bevaEBx+E\nxYvhoousnDfbGT4yU+Xoji9HH+hualVR/PqryIsvmpy33bqJtG4tcu21It9+W+yywmxGF5AjYPeJ\nbEbllp9fLJmCr0hLE7ml4UrJt1Qy3+SQEJPRqqI5mclr7xarxP/DLv+o86UcDo8S+6RnvN0zx2zY\nIPLyy8Wfs1ol5d1VcsUlVrn5ZpGMjNN/VZad4u7YTY7mplbKzURMlojS/P23qTbfqRNMnAjz5sH4\n8aaMXhGF5zpfZxRHqcliepLIU/S0pPnfep6IKQ3YoQMkJHi7N8WkpUFsLNwwNppK7UqUNKxoQkOh\nSxfqR4by/gwL9y64iR71/qTvotGkpJy8xnbuco1edeiQ2ZjVvfvpBGCYsoajngzl+v90YdzEUD7/\nHBo3Pv2ywvPRcXFTiI1NJC5uyhnLQT4zU+Vo9Hb0gY6MVSAoKBD5/HORK68UeeCBct2q6CfxOhyQ\nfnwib9XqJMev6CRSrZrIihUu6rSbLVtmvh/t25tcv65I6usiW7eKRESIvP76ySfOUdKwosrPN9+f\nCy8UeeJBq+RHn0wZ2aGDqWrkCw4dEhk7VqRePZH77hPZt09EzI9xwgRTUeuhh0T273e+CXfk3UZz\nUyvlYrm5IjNmiLRpI3L55SZTf2n1ax101tzEOTmmzdLMm2cKq/pC0LvvPpHISJEPP5Ridf18wJYt\nJun/G294uyf+Ye9ekWf7rZQTmDrN9uBgU43ilVdEDhzwXsesVpOVxWIxCTysVjlxQuS110QaNBAZ\nNkxk27byN+OOvNsajJVypYICkbZtRXr0EFmyxLtB0G4X+cc/TJSpX19k4ECRV18VWbPGO/3asuXs\nHxq8aPNmkeZN8uT9yeUYKlVEVqtkXdRBci0hkhrcXubf/bWcGBQnUrOmyODB5t9/QcHp6keemGFY\nudKM1EHsISGyaMKqUxWr1q51XTO+smasST+UOpfMzOKLUL4gPR2Sk01d1r//hh9+ODPXc0GB+b0S\n7MTpRZvN3PO776BhQ4fXgguzOWVk2AkP91w2p82bodd1eSxpMIzIK+uZc82q7D+Pk5m8fs9py+S3\nQlmyBO4fcpixDT+kzvx3YfBg+Oij0xnKkpPLl0DEZoO1a2HPHvjlF1PTusTfyzXdsG/YwJbgKO5r\nm0zilFBiYpxv8mxcnYHMmaQfOjJWKhD9+adI5coiF19shhKjR5tdqJ9+evbru3cX6dhRpEYNMxsw\nebLIxo0ONeutPMebNom0CD8h2y4bKNKnj8jx425tz1+U5+exfbtZrq1TR2TIEJENby8/NVIttiv9\nq69MIegdO07P0pxrBP3xx6cXfMHsk3j55WLLP+np5qnrOlnl9qar5MvZVp9YnSkrdJpaKSf9+KNL\n1oJ9yvHj5ijIV1+Z32yjR4tMmVL6tYcPiyQlmeNZWVlON+mNureb/7DKLfV+lPR2fUT69zfr7gGg\ncF9BTIzjNY8LueLnceSI+WfTOtwqm6t3kPzgEMlrW6SQxjPPiFx3nVnfrVFDcqLaSlZIZcklSLbX\naiBpf6YUv+HYsSJ33CESFFQssG/ZYk4Gdu4sUreuyIgRIl9/LZKX5/Db9joNxko5qrCqfLNm5uO4\nKhd37Ew9l+3rrZIa3F7sWETCwry74ciFXDXD4MqfR26uyKfvWeXBy1dJw+pWadVKZPhwsyP7119N\ngaTta9fJU/V6Si4m0OYQIgPDh57Z75PVpgqCQ2R3ww5yVTur1K8vcu+9IosX++R2BIc4E4w1Haaq\nuKxWGDbMJC3+5Rdo0MDbPfJ7nqz5m5EBj/ZK4WP7BiwIHD9uFo7r1nV5W5529rOvU85ev7cUrvx5\nhITAbSNCuW1EF14tMEvHP/9sHm+9ZUpAV6tWm5wDH3E719GGDWwkiu8yppHa52+uuiqCo0fhyBE4\nejSUY4eSaVEzldY3teXZoaF06xYYGTqddr5oDUwH9gJ/FnmuNrAY+Av4Dqh5jtd75JOIUg7ZtMlk\nyLr//tM1T1W5eWrNeO9e8+N7aeK5axD7K1eNaD25hm8Gu+8JiNTAKleySmpgFRCJjFwk77wj8skn\nZnn5559F/vrr3KfiXDFN7y24aWT8PjAN+KDIc08A34vIfy0Wy+PAkyefU8o/JCTA2LFwzz3e7klA\nKcx4lJAwpcjOVNcWwDh0CHr2NFUZHxkfCo8kn7MGsT9y1YjWEz+PQqGhEB29g3XrsskilJ/pcqrf\nXbqsZuTIXud8fVEVqphKobJEbCCC4iPjTUCDk//dENh0jte6/2OIUo7yp62Z6pSjR0U6dRJ59NHA\n/hF6a1d6ebmq397YCOhKeHDNuL6I7D0ZafdYLJb65fxMoJRnlTyXq3xedjb06yuMl6fpO2YkFku4\nt7vkNp4c0bqSq/qdkWGn+KwA+GrpwzOIc3k1XLWB65ytT5gw4dR/x8TEEOOOU9tKqYCVkwMDbhYS\nDo7h+rCfsYSVrFEbeFq0iHBos5avcEW/PbkR0BWSFi0i6YMPzAbCLVucukeZMnBZLJYIYIGItD/5\n541AjIjstVgsDYEfRKTNWV4rZWlDKbf54gvo1cvU1lV+JzcXBt4q3LP1cfpXX4bl+++hVq0yv95b\nGcGU80pbM46M9ME140OH4I474Kef4PLL4cYb4cYbsXTogDiYgauswbg5Jhi3O/nnF4FDIvLiyQ1c\ntUWk1A3u2D7HAAAgAElEQVRcGoyVV82cCYmJsHKl76W1VOeVf9hG4m0pXLbjC26tugjLD8scOrrk\nN7/U1RlcnaLSLUTgyy9Nrc6aNU897Uw6zPMGY4vF8iEQA9TFHHFKBL4APgWaAunAIBE5cpbXazBW\n3vH11zByJCQlQevW3u6NclDeIRsZLbvR5GgKlSoHY9mwAVq2dOgew4ZNZM6ccZSc7oyLc+y8rqrg\n/vrLBNuGDct0uTPB+LwT8CJyh4g0FpHKItJMRN4XkcMicr2IXCIiN5wtECvlNT/9BCNGwFdfaSD2\nQ7m5kHBLCuHWVIIpwGK3m0ryDvLrjUDK+7Ztg/h4uOYaWLfOrU1pBi4VeP7+G269FebMgc6dvd2b\ngOKJ9decHLjtNggLi6ZSdFvYtAGiosw5Ygd5fSOQiPklHhxsUlgVfZRxlKW8YMcOeOYZmD8fRo0y\nm7KKTEO7g5ZQVIEnJ8fk6Lv2Wm/3JKB4Yv312DG45RaoXRtmzYKQHFu5Enp4bM04I8PsSSh5ZC4v\nDzp1Ml+LPiwW2LnTde0r1zl40Hz4u/tuGDcO6tRx+BZuWTMuLw3GKpD4685cV/Tb3euvWVlwc78C\noi7cz8tzGzpVirk0btkIZLWavQjffw9Llpgp9LVroWnT8t03Pd3UYe7Vy0yNVq5cvvsp52RnQ/WS\nyxtl50ww1mlqpcrIX1P0uarf7lx/tVqh/435vLg3ns7NKhEUPLPc9yzk8vO6//ynWQLp0gWuvx5m\nz4aOHV1T5eCCC0wQSEiAlBQTkAcPNsdnQkLKf38f4fMfassRiJ3maMouRx9oOkwVIPw1RZ+r+u2u\n93/okMhVV5yQ31sOFPsNvUSys8t1P7dLSxM5dsz97Rw6ZCor9Oghkpjo/vY8xGdSfWZni8yZ45Zb\n40Q6TN9MZ6KUIz76yJRBdDN/3Znrqn5PmhRPZGQiZkMUFK6/TpoU73TfDh6EPrHHeefALVzavgDL\nV19CtWpO388jmjf3TAKZ2rVNNYzvv4f//Mf97XnI2ctDzvBcJ5KSoH17WLDArOH7AA3Gyr99+63Z\nZHHihNubOr0ztyjfTdFXyFX9Lsw7HBc3hdjYROLippRrin7vFhuPdk5m9qHetOkShuWTT3xnjfTY\nMZg61awd+oKzLaAfPOjZfriAVz/UWq1w332mjvlLL8Hcub4z/e/oUNrRBzpNrdxl2zaR+vVFVqzw\nUHM+Mr3mIF/s96/LrJIS3EHyg4LF3ry5yOHDXutLMfn5ItOni4SHiwwaZAon+6pt20Rq1zY1udPS\nvN2bMvPacs+6dSJNm4qMHOn2f284MU2tu6mVf8rJgauvhuHDYcwYjzXrFyn6SuFL/Z4+HT4bt4oF\ntu5UKsg3I5Ply82GKG/66y+zWSosDCZPhiuv9G5/ymLvXnj1VXj7bbOx7MknoUYNb/fqnLyWojQr\nC3791aSudDM92qQqjpEjzf9cc+dqOUQ/ceIEPPww/PgjfDnbRqu7u8GGkwk9kpOdOkfsMrt3Q4cO\n8Oyz5t+Wv/2bysgwgXjpUvjmG7O724f50odDd9BgrCoGux1eeAEeesjnRwHKyMyEwQPzaRu6g//O\na0lYGGArX0IPl9u5s/znhL3t55+hXTvf3wQX4DQYK6V8zk8/wQMD9/JV1cE0i43E8t50b3dJBbIf\nfoAXXzQ7pb20OcsthSKUCgRpaekMGzaR2NhEhg2bSFpaure7FPBE4M034dn+q1ldcAURw7tjeedt\nb3erYvrrL7PPIpAdO2Zmy4YPN199ZZd0GWkGLhXw/DVzlj/L2W9javx6gn//ma8szxP83nTo39/b\n3TLmzoXIyIpVROStt+CLL8yI8fbb/W9N/HxWrjTVlTp3hj//dCqftLfpNLXyD7m5JlWgE7xR03bv\nXti1C2rVMr8XataEoAowD2W3w+cf2Gjzz260ykul0gXBWFathEsv9XbXTHKHxx4zZTU//9wkfahI\nfvgBxo4168kvvxw4H0ZSU01a0jfeMNXafIDmplaBadEiU84sOdmpT/TuTDKQmwubNpkqeX/+ab6u\nWweVc23E1EthTV40O4+EkpVlAnLt2iY416kDjUNtdKiUQnivaDpdF0qzZs4NWHwlz+8PP5hY1y4r\nhQEFqVSSfLBbPJKQ5byOHIEBA0wg+u0384OoaGJjzXv/4ANTGuvOO+H5573dq/KLijK78v38Z6rB\nWPm2XbvM9NPHHzs9tebKmrYisOhTG7++n8LCndGs3RpK8+ZwefQJYhpuIq79elpG/E7ox9OxbLGZ\n6dCda8ivGsqRI3DokHlYvl9C1Gv/pNrBnexYGMV1lVeQExJK165w1VXQtStcdtn5E1L5whR8Sgo8\n/jhs3AjPPQeDbowm6Nq2p48tOVGH2KUOHIAbboBu3UzWJVcUdPBXlSrBiBFmqnr7dm/3xjUsFr8P\nxIBm4FI+LC9PpHt3kWeeKddtXJGBym4X+eILkWs6WGVjFZM5KvviDpK9x2ouuPtukbZtRQYPFrn3\nXpFKlUxqoZAQkVWrzrzhyJEiFsupFET2du3kyF2PyGev7pQHHxS57DKR6tVFenS2ymtDV8rCT6xi\ntZ55G28Wr9i5U2TECJGetX+Vpbf/n+TkFPlLq9W879I6XUbbtm2XuLgJEhMzXuLiJjifMezbb0X+\n/W/zQ1T+Kz1d5LPPvN2LMsGJDFwajJXvGj/eVKzJzy/3rQp/scfGOvaLvaDA/P/foYNIbNu98tuD\n74k9OPjMQFv0F73Val4QEmK+lhaQil4THS2yYIHIc8+J7Nlz6pKs3VaxRZrAv6V6O2lY3Spdu4r8\n5z8iy5aJHD8uEhMzvkQgNo/Y2PHl+Xad08HtVnlz+Eq5JmytrGszWAoaNhJ55x2XtuGLKTwrhKws\nkw7UE1Wpyio7W2TCBJE6dcz/I35Ag7EKHJs3izRqJJKZ6ZXmCwpM9bp27UT6RafJtj4Pir12bfNL\n4XyBVqRsI8PzXbNypUhh4AfJj4qWtKFPyPuDF8p1naxSo4ZIw4ZbBU5IDazShZVSA6vLR8YnTogs\nXy6SkCByXSerbLBEST5BYg+qZEacNpvD9zzfqNdfy1X6vfR0kX79ROrVE3n8cfNnb8nOFpkxQ6RZ\nM5MnfLv/fBDTYKwCy4EDHm8yP19k7lyRqCiRwdEpsjN2uNjr1DG/mHbvNhe5YAq2TIqOntu1M6Pn\nhAQzdV+9uuTGXC/vvLNXwsMWyhraywmCZQ3tpUG1lfLQQ4dl3jyRv/5yfGLBbhdJSRF55RWRvn1F\nwsLMtPnjj4v88tpKsQcFnXsK/jzKMur1xohfFbF5s8iYMaYQxa23iqxZ4/k+jBgh0ru3SFKS59su\nJ2eCsR5tUuqkAwfgH7faaHgghSGT2nL9tJuw3HADPPCAOaPkDWdLGXn8OGzeDO3bk/nZfC68/XZC\nxE6epRLfjpxJfmYuq45GsXhnG7buD6N1a7jiEhuRx1PYERbNsUqhFBRw6pGfb74GH7eR8+ufVAku\noEPXalw89Aquuw7q1SvSn27lyyldlqNmTh9HS0kxP8iYGIf6pM7CZoNZs0yu66uu8mzbdrvfngd0\n5miTjoyVEpH160WiI6yScWEHsybcoYPI0aPe7lbZlFyj/u03kWHDzHC2WjUpaBwuhy+PlaxajSU/\nKFj2N+kgM1+3yvvvi3zwgch3T6+WHVcOlMxLb5TcKjXEDmKvXFnk1VfP3l45ZgbKMup1as34999F\nGjQwUxvKMz7/3PzPk5fn+GuPHDHr0wkJru+Xl+HEyFiPNqkK78svTaGeWQ+k0Pi5VDNM3LDBPLxd\n1q8sQkPNCLXoCHrWLPN3djtBO3ZQa/58eGw52Auot3cDd16eevq97QqHNkNg61b49+JTrztrUojQ\n0HJ9X8py1KxFiwiWLBlNQsKUIpV9znFca/VquOkm+N//zBla5X4i8MknpixhZiZER5vkLpdeCvfc\nc+ao9uBBuO022LPHPHJyoFcvuOsu7/Tfx+g0tfINR4/C339Dp04ea1LEVMz765WFPPLRlVx2ZUi5\np2B9Vlmml10wBV0WLq9nu3KlSegxYwb06ePi3qoysVpNtps1a8yHuldfPfOavDxTt7pRI2jY0JwN\nDrS0nCdp1Sbln0RMUfcGDWDaNI80mZ0Nd8cXcMPKRO7kA4IXfWNKz/laWT9XKst789D7d1k926ws\naNMG3n3XjLKU8gEajJV/evttk1f255+hShW3N7djB9zZ9yCvHYqj7UUnqPTpx1C/vtvbVW5y8CDU\nrevtXih1ipZQVP5n/Xr4979NuksPBOIVK+Cey//gy4wraDe0HZWWLtFA7O80EKsAoBu4lPdkZZnp\n6SlToHVrtzf38bs2Zv4rhTd6LaDmwP+a/LxKKeUDNBgr79mwAXr2NNVj3GzhJzai7+/G16QStKkt\n9H7S7W0qN8jPh2D9taUCj05TK+/p3NnsunTzjspVq+C1e1NoQypBhceWUlPd2qZygx07zCa7nTu9\n3ROlXE6DsQpoGzeaY6ePTI8mqG1bCAnxjbJ+yjGZmdCjB9x7LzRt6u3eKOVyuptaBaxdu+CFjh8R\n+1RXBo6NCOxjS4Fs3z649lr4xz/giSe83RulzkuPNil10uHD8Gz7j0m0jiX0tyS4+GJvd0k548AB\nuO46uPVWmDDB271Rqkz0aJPybf/9r0mf52bHj8Pz13xD4qGHqLF8oQZif5aaarJrJZ6jOIRSAUBH\nxsozli2DuDiTx7ZJE7c1k58PE2KSeOy326nxw9cEdb3SbW25W2GWqowMO+Hh5chSpZTyKJ2mVr4p\nMxOuuAI++ACuv95tzYjA43G7eGreZVT/+mNCboh1W1vu5vL8zUopj9FpauV78vJMYo/773drIMZm\n4927V7F6Y02Cf1np14EYICFhRpFADFCdrVsnkpAww4u9Ukq5iwZj5V5PPw01apiUl+5is3Egqhvx\n73dnaW43akQ2cF9bHpKRYad4iUGA6mRm2r3RHc/YsQO++87bvVDKKzQYK/e66y6YPfvM2qYulPJR\nCmG7Ugkhn5DNgZHQ43TN36KK1/wNKNu3Q0wMbNrk7Z4o5RW6Zqz82sGD0K2jjVXB3aiZETh1iL2x\nZuy1DWPbtpnjS+PGwahR7m9PKTfTDVyqQrHb4a12b2Dr3IPHXwsPuIQeLqv5W8a2vLJhbPNmk1nr\nqafgvvvc145SHqTBWFUon965gO4fP0CdLb8S0rSht7vj14YNm8icOeMovk6dTVzcFGbPdtMZ3/x8\niI42I+KRI93ThlJe4Eww1vInynUKCuCzz0xpQjcXf/hl9t/Ezr4b+xdfaSB2Aa9sGAsOhuXLtZ60\nUugGLuVKjz0Gb7xhRjxutHeLjZojbmHv6Geof1MXt7ZVUXhtw5gGYqUAnaZWrvLGGzBtGqxcCXXq\nuK2ZgnxhReNBVG9Siyv+eMel967IGa80yYhSrqNrxso7vv4a7rkHfvoJWrZ06a1LBsjatR+m5vJk\nJq66gUrVKru0nYoejNy6YWz7dvjjD1PwQakAp8FYed6aNXDDDbBgAXRx7ZTxmQEyh0qVTrBypZXO\nnV1b09YrG5gqiqQkGDrUJH7Ro0uqAtB0mMrz6tc3ST1cHIihtJSQVSgouIDXXnvP5W1VyIxX7iZi\nli8GD4ZZszQQK3UOuptalU94uHm4QWGArIGNaFJIIZosQt0SIE9vYCo+Mg7YjFfuduKECb6rVpl9\nBJGR3u6RUj5Nf9MonxUeHkQNdvMznfmR7iTTjRrscUuAnDQpnsjIRE7vKDZrxpMmxbu8rQphzx7I\nzjbBWAOxUuela8bKZ6WlpfNVu/sZnb2QIOAEIcSF38bk5OeLbSxy1S7ocm9gOnzYbFL64w/4/Xd4\n7bXSj+48+ywcOwbVqkGtWtC+PVx6qSmooZTye7qBS7mX1QrvvQcPP+z2pB4Ah3/agP2abmRXvYBG\nOfvIrHkhsnwpzdu1PXWNT+yC/te/TLKT/fuhY0e47DK4/HIYMADCws68/p13YN8+E5APHjSb4FJS\nID0d6tXzTJ9dScQj/x6U8hcajJX7HDwIN95oAs2bb7q1ChOAHDvOjkadWX3lGAZ/Nuisead9Yhd0\n4dnqiy+GSpWcG6nn5UFIyJnP5+aaylfXXAOxsdCqle8EPqsVXnjBHFv68ENv90Ypn6HpMJV77N5t\nji/deCO8+KJHgsFfN45hW1A7Bnx1F1SxnHW3tkd2Qdvt8NVXYLPB8OFn/v1VV536z9JG6qtXl2Gk\nXlogLmy7Rw9zPOi550x2s5gY6NMHhg0rx5sqh7w8M7p/+mnzb2LyZO/0Q6kAohu41Lmlp0P37jBk\niMcCcVryLk789BstFr1F5Srnbs+taRzz8syRnHbt4JlnoG7d877kzONY1dm6dSIJCTOc60OVKjBi\nBMycaX4WK1aY4Lxnj3P3K69HH4VGjeDzz2HRInj/fbftpleqItGRsTq3UaNg9Gh46CGPNJeXB0PG\nNSFuyq88dOX5A+qkSfGsXp14xprxpEmjne+ECLz1Fvz3v9C8Obz8MvTsWaYPIm4dqVssJsPZubKc\nzZ1rPjR17AgdOpz+6miK0rOtA992GzzyCDRp4tj9lFLnpMFYndtnn8EFF3isuUmToHZtGP1w2Ua2\nLVpEsGTJaBISphTZBV3OzVsWC2RkwJw5xaagy8Lr55UHDDBHidatM4/PPoM//4SxY2HChDOv//xz\nmD//dOAVgbVr4R//MKUNS+ra1a3dV6qiKtcGLovF8ghwN2AH1gMjRCS3xDW6gauCcfao0U8/wcCB\nZnNxo0bu76c7+MTu7pLsdsjJMUepSioM2mACMUCbNmajXrB+VlfKGR7dTW2xWBoDK4DWIpJrsVg+\nBr4RkQ9KXKfBuAJxNhhZrWZG9eWX4eabPdZdt3BrwQWllM/zRjBeBXQEbMDnwKsi8n2J6zQY+4Oj\nR80O2UcfLdcmLaeOGu3axbuDFvNL9F28/bbTTTtn3jyTdKNVKw83rJQKVB4tFCEimcBUYAeQARwp\nGYiVn/jxR7PJZ9s2s4OqHBzewHT4MIev6cfxv9J5+eVyNe2Y3FyTvOSxx0zyDaWU8iKnF4UsFkst\n4GYgAjgKzLNYLHeIyBmn/ycU2TgSExNDTEyMs80qVzpxAsaPN8d33nkH+vYt9y0d2sBks1HQqg01\nD+zj7ouEavZxQOiZ17najh0waBA0aGDSVtau7f42lVIBKykpiaSkpHLdozzT1LcBvUTknpN/Hg5c\nKSKjSlyn09S+KDPTJI5o3twE4gsvdMltHVkzLnj5NYLGPowFTNKL5cvdUoqxmIULIT7e7BQeN853\nslkppQKGp9eMOwPTgU7ACeB94FcReaPEdRqMfVFeHnz5pdm+7OKAVKYNTLt3Y72oIwdzqtBEdpWa\nd9otZs2CiAiTyEQppdzA47mpLRZLIjAEyAPWACNFJK/ENRqM1RkWv7WaeQ/8xly5nbakkUpzGkRO\n8e4RIKWUcgEtFKFKt29f6aX8vMRmg8aND5GVVQ2oUuRvPFzgQSml3MCju6mVjxMx66MxMWZjlg99\nIBo9GmrUSKN4IAaXF3hQSik/ocE40BQWN2jfHp54Au65x5T485GNSnPnwqpV0L37d7itwEOhmTPh\n559ddz+llHITnaYONHFxsHevOT9bxuIGHlFQwPYdQXTqbGHRIqhTx81pI2fNgiefhB9+MHWGlVLK\nQ3TNuCI5W1WdY8dKz0HsZfbHn+T/Pq7D8VH/OlV/wG1pIz/+2FQWWrrU5FlWSikP0mAcyLKzTSWF\nH36AZcugdWszDesPFi7k6JB7GXnpH3y87EKC3Lk4Mn8+PPAALFli6hArpZSHaTD2Z2cb6W7YYKae\nN2+GSy+F2Fi47jqTHKNKyQ1QPmjnTnI7dmIQn/J/Kd3cW41pzx64/HL4+mvzvVJKKS/QYOwPDh6E\nV1+F/fvhwAHzddcuqFHD1JEtKTsbNm0yI+HqJXM++zCbDdauJf+RfzF16wDaznqCfv080O7Ro1Cz\npgcaUkqp0jkTjLVgqSscO2YC5saNZiS7caMJCkuXnnltpUoQFATR0SYFZb16EB4OzZqVfu/q1c1o\nz5/YbNCtG7J+PblB1Tkw4jvPBGLQQKyU8ks6Mi6vY8dMQI2MhKgos2EoKso8oqO93TvvMGeXID+f\nPEIoSFpOlWvdnHNaKaV8hE5Tu9OuXSbolrZOW1BgRrzKsNnI6dQNy18bsLeKoupvyRDqgWpMSinl\nAzQDl6tZrfD++2bDVIcOkJJS+nUaiIs5UhBK1/xkliYsd28gnjjR7DBXSik/p8G4NL/+CoMHQ9Om\n8NVX8OCDkJEBV1zh7Z75PLvdbP7u3jeUPk93cV8gnjMHPvjAbGxTSik/pxu4SnPsmDlC9OabULeu\nt3vjH+x2mDKFZw4+SFZWdaZMcWNbv/8OY8aY89b681FKBQBdM1au8cQTHFiwiiutS1j1+wXuKxK1\ndy906gSvvAK33uqmRpRSynl6tEl5x4wZ5H74Kdcc+5mPFroxEIvAHXfAiBEaiJVSAaXijozz8mDq\nVMjKgmee8XZv/FdyMvZbB9KvehKDJkQRH+/m9lJSzLExt+bUVEop5+nRprJaswbuugvq14f//Q+a\nN/d2j/zTjh1I584kRszkUKdevP66tzuklFLep8H4fPLyIDER3n0XJk+GO+/0nRKDfiItLZ1nn/gf\nlTfvIb9lPdrlX8UnhwawdCmEhHi7d0op5X26Znw+U6fCb7/Bn39Cw4be7o3fSUtLZ0CPqcxMSyKK\njWxY04aYoNEsXrWTkJCm3u6eUkr5rYq18Pbww7BggQZiJyUkzKBa2gCi2MgF5NOGTbS2b+O1195z\nT4N2O/z4o3vurZRSPqRiBeOqVaFyZW/3wmelpaUzbNhEYmMTGTZsImlp6cX+PiPDTgqd2EBbThDC\nRqJI5TIyM+3u6dALL8D48SYoK6VUAKtY09TqrNLS0unZcxpbt04EqgPZrF6dyJIlo2nRIgKOH2fc\n3iUM4X66kUxbUkmlLVkE0bixGz7TrVgBr71mlhV057RSKsAF7m85mw2OH/d2L/xGQsKMIoEYoDpb\nt04kIWGGqancty9XtWyA/YIssqjCz3QhiyAiIxOZNCnetZ05eNCcJ54+HZo0ce29lVLKBwXmyPjo\nUbjxRvMLfdQob/fGL2Rk2DkdiAtV58iO49C7N3ktLmbAtne4eeAxLJbn2L3bTuPGQUyadHLk7Coi\nEB8Pt98Offu67r5KKeXDAi8YHz4MvXrBlVeaAg+qTMLDg4BsigbkmmTwyobZ5N/Snxu3vkFkqyDe\nfTeUoKBE93Vk925TpvL5593XhlJK+ZjAOmd84AD07Ak9ephzxHqGuMxKWzN+pXZPhvdrzdA907mw\nvoWZM7VapFJKnU/FTvpx8CDExMBNN5n0lhqIHVaY0KPKlj3kXNSQfyWO5NF/taBadQsffgjBgTeP\nopRSLlexg/GxY/DllzBkiAZiZ9ls0K0bpKZij2rLsGbJHA8O5ZNPNLuWUkqVlTPBOHB2U1erBkOH\naiAuj5QUSE2F/HwKUjbQ8GAqH32kgVgppdwtcIKxct7hwzB6NISHI1FtyQsKYUe1KJ77sq37c6Ss\nX28SeyilVAWmwbii++EH6NABLBb22i/k5jrJPHr5chptSabKhaHubTs7GwYNgosucm87Sinl4/x3\nzTgry5xJDXVzwAhUJ05AQgLMmQPTp/Nj1d7ccQeMGAETJnhos9aIEeZnOGOGBxpTSinPqDhVm0RM\nPeLWreHpp73dG/9hs5l14chI6N0bmjXD/sdaXph+Ia+9ZmJi794e6susWbBqlUl3qZRSFZx/BuOX\nXoK0NPjgA2/3xO3S0tJJSJhBRoad8PAgJk2Kdy7jVZGd0rRtC5Mnc6Dj9Qy/04LNZmKixzJP/vUX\njB0LS5dCjRoealQppXyX/wXjpCST0OPnn02mpgB23uINjiiyU5oNG1i3LZT+d1sYOtQcy/bojun6\n9WH2bGjf3oONKqWU7/KvNeNdu6BzZzMivv5619zThw0bNpE5c8ZRPGd0NnFxU5g9uwwpKfftM6PP\noUNPjYxlwwb214uiS14yr70fSr9+7uq9UkpVTIG/Zjx3Ljz00BmB2GVTuT7mbMUbStYPLsycVXnz\nHk5c3JAJI2Jo8unHMH++SYIyeDCEhrJ3XjKT41NZk9uWpHmhNGvmsbeilFLqHPwrGP/rX2bzVhEu\nncr1MaUVb4DsYvWD09LSGdBjKjPTkmjLBo6vqULOvFc5POZ+am/eDPXqsXOnmdmfPTuUe+/twsKn\n4YILPP1ulFJKnY3/nTMukWHrnHV4nZCWls6wYROJjU1k2LCJpKWll6u75TFpUjyRkYmYgAyQfUb9\n4ISEGVRLG0AUGwmhgKrkMNA+n9F7Q9lmrce995pjxFWqwIYN8MILXgjEeXmQm+vhRpVSyn/418i4\nFGWdyi2Lso6yXTotXnjcKDr6jDPTLVpEsGTJaJ594lkabNxKk+q53NGmDjV/WAot7gLM+0+hExto\nSxs2sJEo1nIVm79vwKJF8MADsHkz1K3rXPdc4qmnzNfJk73YCaWU8l1+H4zLMpVbVmcfZZ/eMOXQ\ntPg5Au2pv7/6ati40Rw3Sk4uft1339Fi7Fje/esvKCiAsDC48VG44opi7z+JILqRTFtSSaUtWdQg\nsuEhkpKgVi1O9dsr6+rffAMffwxr1ri/LaWU8lci4taHacJJCxaIrFt3zku2bdsukZGPCmSJWVDO\nksjIR2Xbtu0ONxcTM/7kPYo/YmPHn7omLm5CkbbkVJtxcROK38xqFenQQSQ42Hy1WkV+/11k4ECR\nbt1ELrlEJDT09E1CQkRWrSp+j337RGbONPco5ZqcHJG5c3dL7dqLBfJP3ipHWrR4stj7d+X3yCE7\nd4o0aCCyfLl721FKKR9yMu45FCt9d2ScmWmybH333TkvK5zKTUiYQmamncaNg5g0ybnNW2UZZRdO\ni1Bk0pcAABPQSURBVNfARjQppBBNNY4Rvv5vsyC7dq1JZHH33cXO9ZKaCk2bmp3N9eubR9WqcPPN\nZmQcFWVGx0VdeCHccotJcrJhA9ImihR7WxZPhe+/h59+gtatGzJ4cFU2b55Dfv52mjQRJk36Z7H3\nX5YRv8vl50NcnClA0a2be9pQSqkA4ZvBWAT++U+4/3649NLzXt6iRYRLgsqkSfGsXp1YbArabJga\nfeqa8PAgarCHZHrTllQAsqhO5t7asL8h9O1rppGbNDHBdcOG04E2NBRuv714oytWnM6KVWSKOifH\nJBnbujWUHUOT2f19KrPXtCUkPpTrr4d77jFppevUAagJ3HnW9+XKdfUymzXLfNh44gn3taGUUgHC\nN4PxrFmwcyd89plHmy3LKHvSpHhOJI0lKiOVEPLJJZh7Gt/I5BUvQMnReHJyqYEWzOeN/fshPT2U\nrdu7sHUpbNsGW7eax/790KyZSSMdGRnKZUO78ON0nDob7Mp19TKLjzezAJUqua8NpZQKEL6XgSsj\nw4yGFy+Gjh3d17Gy+vtvePdduO8+aNkSgO3rU7F070Hjo/vJrHkhsnwpzdu1PeOlBw/Cli2wffvp\nR3r66a9VqkBEhLmtCbqnH02bui6OlbbpLDIyMM5iK6WUr3EmA5fvBePvvzfrruPGua9T57N/P7z5\npunL33/DP/4BjzwCjRqdvsZmKzbqFTGBd8UKs5a7YgXs3g0XXwzNm5/5iIjwbPXHwt3Up0f8gZGl\nTCmlfE1gBGNvW7QI+vc3R4maNYPffy/1kG5urjmtUxh4f/rJJNO45hrzuPpqaNfu3KPbQE3jqZRS\nFVng56b2hIIC81XE7OgukTFj1y547TUzc920qQm6t90Gr7zi2HpuQKXxFIGJE+Hee6FxY2/3Riml\n/I7/pcN0t+7dzdRzSEix40Zr1sCwYabqX26uGTCvW2dms++4w/GNVa5O4+lVU6fCggWnM4wopZRy\niI6MSwoNPbUL2t6mLYuSQ5k6Ff76yxSMev1118Qcrxw3codFi8w56J9/hmrVvN0bpZTyS94fGe/c\nCW+/7fl2Fy825QVLWc/OCQllemoX2l0VylNPmVM627bBY4+5bvB3+rhRUW4+buRqf/8Nd94Jn3xi\n5uyVUko5xbu/+UVMJYM9ezzTns0GS5ea7FgjR5oMXyWqQC1ebI4WzZtn1obXrIHhw11f6agsFZl8\nWna2yR727LNmx5pSSimneXc39fz58J//mIhXubJb+4HNBp06mfnm2rXNgm+R0Vx+PowfDzNnmpwj\n113n3u6Anx83EoFly6BHD2/3RCmlfIp/HW2yWs0GqQ8/NJum3O2rr8xIDszmrOXLoUsXwMyUDx0K\n1aubQFy/vvu7o5RSKjA5E4y9N02dkAC9enkmEANcey1cdNEZu6S/+cYMmPv2hYULNRArpZTyPO/s\nps7LM/kgp0/3XJs1a8Iff5zKmpVbOZSnxpm9R/Pm6bKnUkop76mQGbi2bzcbqevVgxkzzFdXC6js\nWvn5Ziv5ww+bPJ5KKaXOyuMZuCwWS03gXSAasAN3icjP5bmny4icsVMa4PPPTXXGJ54w6aZLuaTc\nAiq7Vn6+yc29f7/O4SullJuUd834VeBbEWkDdAA2lr9LLvDnn9CzpwkkRbz9thncLVgAY8e6JxBD\nAGXXKgzEBw7Al1+a+sRKKaVczumRscViCQO6iUg8gIjkA1YX9ct569ZB797w6qsQfPrtffSRSZ/8\n449mH5c7BUR2raKB+IsvNBArpZQblWdk3AI4YLFY3rdYLH9YLJa3LRbL2X9j//IL7N1bjubOw2Yz\nh4R79oRp02DQoFN/9e23MGpUAR07/h/33JPIsGETSUtLd1tXAiK71sKFGoiVUspDnN7AZbFYLgdW\nA11F5DeLxfIKcFREEktcJ4lPPGEqKtxyCzHx8cTExJS748UUTegREQHr158qFpycDAMGFFC16ptk\nZNxF4RpuZKT71nBLWzN2Z3tuk59fbHZBKaXUmZKSkkhKSjr154kTJ3ou6YfFYmkArBKRlif/fA3w\nuIj0L3GdyJgxcOQIvP++U22d16pV5myS3V4soccff5gZ63btZrFs2a0UnzrOJi5uCrNnJ57truXi\n19m1lFJKOc2ju6lFZK/FYtlpsVhaicjfQA9gQ6kXf/ihOd/rLtHR0K4dbNhwKqHHpk0mkcdbb8G0\naVvw9BpuixYRbgv0SimlAkt5FzEfAuZYLJa1mN3Uz5V61X//6/Rh3rS0dIYNm0hs7DnWegvLHi5f\nDsnJpB8K5YYb4Pnn4dZbA2QN110OHzaVMDb6xkZ4pZSqiDyT9MNud+ockTNrr3v3Qrdu8OCD5hiT\ns/epEJKSTAnEW26BF17QjVpKKeUC/lUoogyGDZvInDnjKOta75EjEBMDAwbAhAnF/07XcIvIzTUl\nqmbNMilJe/f2do+UUipgeDwDl7ud7bzu3l15ZkT3/PMQHg7AsWNmjfjaayGxlKVaXcM9SQRuuAHC\nwmDtWrjwQm/3SCmlKjyfDsan13qLj4wf2bsMqtWCRo1OPTtmjClP/PLL7susFRAsFrOr7ZJL9Bul\nlFI+wqenqUtb6x3VYAgvhawlZN06qFMHMFWXnnjCFGUKC3Nd35VSSilHBdyaMRRf672i+kGeXzmX\nSkuWwGWXAaYSY6dO8PXX0Lmzq3rt57KzYc4cWLrU5AHVEbBSSnlMQAbjU3Jz4dJLTYWHu+8GTIKo\nmBi46SZT4a/C27rVZDqbOdMkQRk1Cnr00GCslFIeFLjB2GaDlBQ4fhyuu+7U04mJsHIlfPcdBFX0\nI8MjR/5/e/cfY2V153H8/UVcQRgXSxQXEYbFhE2pLVpW8Udd2LS7TTe1drf1R0tWQmvRsGKqIeCK\nQUvcbLf+iLG11hW2zUaxq11/NNafayfjVEFbhUWEokasMooFtlQw2tE5+8e5IrAMy8zc53nuZd6v\nZMJwuXPPN3nyzIdznvMjn6w0axZceCG0tlZdkSQNSAdmGL/1Vl44vGYNTJqUN/doaaG9Hc4+Oz8n\n3mUe14Hp7bdh3Tr41a9gypQ8QrCnDRvgqKNgyJDSy5MkfagvYdz4/cnnnstB/N57ebvLNWvYuhVm\nzMhLZJsyiLu6YOPG/D+J3/xm7++5/noYPx6GD4fDD89LuTo64J139v7+1laDWJKaVPP0jGv7Tqf2\nx/m7mS2MG5fzqmksWQLXXZe3CNu2LW8POmoUfPOb+dzgPXV25mH5I4/MgexzX0lqCgfeMPWOHTBs\nWA7k2jD1zbe18IMfwPLlcMgh9a21T1KCF17IBa1Ykdfvzp37f9/36qt5i7BRo2DkSDjooPJrlSQV\n7sAK487OvFZp1aocXuQR6+nT82jtxIl1LrS3Vq7Mi5ufeiovbj7pJJg6Nc9e/vjHKy5OklSVAyeM\nU8obTH/iE/CtbwF5xPbEE/Oo7qxZBRTaW6+/nnvCJ53UpA+uJUlFOHDCeNkyuPrqPHu4NhY9Zw5s\n2ZL/qdTHp11deUh5wK+dkiTtj6Y7KOKD3bU2buzm6KNrJykNG5q7vz/96c4gvv9++NnP4NlnSwzi\nzZvhllvyJho//jGcempJDUuSBprKwnhv+04vX76IZydvoOW88/Iel+Q5XHPmwK23wogRJRS2ZQtc\nfnkO4C9+Me+zOXlyCQ1LkgaqysZer7jih7sEMcAwXnrpKha+f8xuhxEvXpw7pZ/+dAlFvfhiDt4h\nQ2D9eli61CCWJBWusp5xT2cVr952GAwdCuTZ00uWwOrVJRXV2pp7xKecUlKDkiRV2DP+8KziXe1g\n9OhcUnd33mL5qqvyLo+lGDzYIJYkla6yMF68eCYTJixiOG8wlScZzhtMmLCIxYtnAvngoXffhdmz\nq6pQkqRyVLq0acPqNQw6/S/5k99tpnPEEaT2/6L1uEls3pzPhHjggZ3HFtfXpk0wbx5cey0ccUQB\nDUiSBqqmOyiidfvvGfv7zRxMN+N2bKV1x1sAzJ8P55xTUBDff3+elDV2bEnTsyVJ2rdK1xnT1ZU3\n0zjoIPjoR2HSJDo68vnEzz9fQHvf+x58+9t5ktbppxfQgCRJvVfdMHV3N5x8cj6x6IQTYNIkuoa0\ncPzxsGgRfPnLdS7k+9/PQfzzn+ejCSVJKkBz7cC1dGnuEV9wwc6tJq//FzjmGPjSlwpob9AgeOwx\ng1iS1HCq6Rlv2ZKHpR96aOemGhs2wCc/mQ9BmjCh0JIkSSpM8xwU0dUF7e35uMGaM87IpzItXFho\nOZIkFap5wngP996bZ1CvWrXzbAhJkppS0y1tAti+HS66KM+vqlsQ33FHQdOxJUmqv8rD+Oqr8yqj\n6dPr9IG33QaXXFLyoceSJPVdpcPUL78MU6bkgyBGj65DY8uWwaWXwqOP5glikiSVrLGHqS+7DH75\ny91eWrAALr64TkH84IO5R/zwwwaxJKmplNMzfuKJvHh47Vo47DAAfvELOPdcWLcODj20n41s2wYT\nJ8Kdd8KnPtX/oiVJ6qPGnU09eXI+mOErXwHy5ltTp8LcuTBjRp0aevFFOPbYOn2YJEl907jD1CNG\n5G5wze23Q0o7s7k+DGJJUpMqp2e8YkXe0QN4++08orxsGZx2WqFNS5JUusbtGX/jG/BWPh7xmmvy\n+RAGsSRJWTk944MPhvZ2OsdO5bjj8qTqfp3XsH17nvk1ZUrd6pQkqR4a99Sm2lnFl8+F88/vZxCn\nlD+kpcUwliQdEMoJ48cf55kXWnjgAVi/vp+fddNNeYnUk0/WpTRJkqpWyjB1d3di2rQ8e3r27H58\n2IoV8PnPwxNPOHtaktSQGnYC1913w9at8LWv9eNDtmyBs86CW24xiCVJB5RShqnnzYObb4bB/Wlt\n7dqc5meeWbe6JElqBKUMU48e/Ws6Og5h/PhxhbYlSVLVGnaYurNzDJ/5zI28/PIrZTQnSVJTKenU\npkN56aWruOKKH5bTnCRJTaS8IxQZRmdnd3nNSZLUJEoM4x2MHt2L5jo74cYbiytHkqQGUVIY72DC\nhEUsXjxz/96eEsyZA2++WWRRkiQ1hFLC+KtfvYZHHrlo/2dT/+Qnee/phQuLLUySpAZQzkERvWlj\n61b42MfgrrvglFOKK0ySpAL0ZWlT44XxrFkwbJjPiyVJTalxT23aX+++m79uuKHqSiRJKk3j9Ywl\nSWpiDbsDlyRJ6plhLElSxQxjSZIqVn0Y33MPbNxYdRWSJFWm32EcEYMi4pmIuK/XP/zaa/D1r8Mf\n/tDfMiRJalr16BlfDDzfp5+cPx8uvBDGj69DGZIkNad+hXFEjAE+B9za6x/u6ID2dliwoD8lSJLU\n9PrbM74emAf0biHx++/D3Lnwne/k3bYkSRrA+rwDV0T8DbAppbQyIqYBPS5wvvLKK3d+P23aNKYN\nHQof+QicfXZfm5ckqSG0tbXR1tbWr8/o8w5cEfFPwAzgPWAo0AL8Z0rp7/d439534OruhkHVT+aW\nJKmeKjsoIiL+Arg0pXTGXv7N7TAlSQOG22FKktSEPChCkqQ6auye8datpTUlSVIzKSeM33kHpkyB\nlStLaU6SpGZSThhfey1Mnpy/JEnSbsp5ZjxyJDz9tNteSpIOeI37zPiCCwxiSZJ6UE7PePt2t72U\nJA0IlW36sc8GXNokSRpAGneYWpIk9cgwliSpYoaxJEkVM4wlSaqYYSxJUsUMY0mSKmYYS5JUMcNY\nkqSKGcaSJFXMMJYkqWKGsSRJFTOMJUmqmGEsSVLFDGNJkipmGEuSVDHDWJKkihnGkiRVzDCWJKli\nhrH2qa2treoS1Edeu+bm9RtYDGPtk78QmpfXrrl5/QYWw1iSpIoZxpIkVSxSSsU2EFFsA5IkNZiU\nUvTm/YWHsSRJ2jeHqSVJqphhLElSxQoL44j4bESsi4j1ETG/qHZUjIjYEBGrIuLZiHiq6nq0bxGx\nJCI2RcR/7/La4RHxcET8OiIeiog/rrJG9ayH67coIl6LiGdqX5+tskbtXUSMiYjHImJNRKyOiLm1\n13t1/xUSxhExCPgu8NfAJODciPizItpSYbqBaSml41NKJ1ZdjP5f/0a+33a1AHg0pTQReAy4rPSq\ntL/2dv0ArkspnVD7erDsorRf3gMuSSlNAk4G5tTyrlf3X1E94xOBF1JKr6SUuoA7gC8U1JaKEfgY\no2mklDqA/9nj5S8AP6p9/yPgzFKL0n7r4fpBvg/VwFJKb6SUVta+3w6sBcbQy/uvqF+2RwOv7vL3\n12qvqXkk4JGIeDoizq+6GPXJkSmlTZB/YQBHVlyPeu8fImJlRNzqY4bGFxGtwGRgOTCqN/efPR/1\n5NSU0gnA58jDLqdVXZD6zXWMzeUm4E9TSpOBN4DrKq5H+xARw4G7gItrPeQ977d93n9FhfFGYOwu\nfx9Te01NIqX0eu3P3wJ3kx89qLlsiohRABFxFPBmxfWoF1JKv00fbgTxr8CfV1mPehYRg8lB/O8p\npXtrL/fq/isqjJ8Gjo2IcRHxR8A5wH0FtaU6i4hDa//LIyKGAX8FPFdtVdoPwe7PGO8DZta+Pw+4\nd88fUEPZ7frVfoF/4G/xHmxkS4HnU0o37PJar+6/wnbgqk3Dv4Ec+EtSSv9cSEOqu4gYT+4NJ2Aw\ncJvXr7FFxO3ANGAksAlYBNwD3AkcA7wCnJVS+l1VNapnPVy/6eTnj93ABmD2B88g1Tgi4lSgHVhN\n/p2ZgH8EngL+g/28/9wOU5KkijmBS5KkihnGkiRVzDCWJKlihrEkSRUzjCVJqphhLElSxQxjSZIq\nZhhLklSx/wVPOfgaAIZ8VAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1c2072ef198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot to compare the true relationship to OLS predictions\n",
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "\n",
    "ax.plot(x, y, 'o', label=\"data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "ax.plot(x, res.fittedvalues, 'r--.', label=\"OLS\")\n",
    "ax.plot(x, iv_u, 'r--')\n",
    "ax.plot(x, iv_l, 'r--')\n",
    "ax.legend(loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
